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ABSTRACT 

The process of gravitational accretion of initially homogeneous gas onto a solid ball 
is studied with methods of fluid dynamics. The fluid partial differential equations for 
polytropic flow can be solved exactly in an early stage, but this solution soon becomes 
discontinuous and gives rise to a shock wave. Afterwards, there is a crossover between 
two intermediate asymptotic self-similar regimes, where the shock wave propagates 
outwards according to two similarity laws, initially accelerating, then decelerating 
(and eventually vanishing). Lastly, we study the final static state. Our purpose is to 
attain a global picture of the process. 
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1 INTRODUCTION 

The problem of spherical infall of gas onto a solid object has 
application in several astrophysical situations. For example, 
it can represent the infall of gas onto a neutron star from 
the surrounding cloud. Alternatively, it can be a simplified 
model for the formation of planetary atmospheres as a result 
of gas accretion onto the solid (rock) surface of a previously 
formed spherical core. 

Aspects of this fluid dynamical problem have been 
treated before in the literature. A simplified formulation 
with constant gravity was considered by Bisnovatyi-Kogan, 
Zel'dovich & Nadezhin (1972) and a self-similar solution 
with a shock wave found. Self-similar solutions with the 
(variable) gravity due to a central point-like mass, appli- 
cable to either ordinary or neutron stars, were studied by 
Cheng (1977). Similarity methods favour power-law distri- 
butions. This property was amply used in the study of su- 
pernova explosion and later fallback by Chevalier (1989). An 
extensive study of self-similar spherical accretion with an 
initial power-law radial distribution of gas was provided by 
Kazhdan & Murzina (1994). They concluded that a realistic 
solution should be an interpolation between a self-similar 
solution with zero mass flux near the origin (of which they 
derived general features) and a solution with constant grav- 
ity and the correct boundary conditions at the solid surface. 

Chevalier (1989) and Kazhdan & Murzina (1994) as- 
sume that the incident gas flow is cold (T = 0) and, respec- 
tively, that its accretion rate decays with time with some 
exponent or that its density decays with radius with expo- 
nent < 3 (both conditions are related). A cold inflow is not 
a good approximation for large distances where the gas ther- 
mal energy is non-negligible with respect to its gravitational 



energy. On the other hand, the condition lj < 3 implies that 
the mass of gas diverges at large radius, making the neglect 
of its self gravity questionable. 

One may try to amend these problems by endowing 
the incident gas with pressure. For polytropic flow, similar- 
ity demands that the initial pressure be also a power law 
with radius and, furthermore, determines the exponents (in 
terms of the polytropic exponent). This is the case studied 
by Cheng (1977) and it may seem somewhat unnatural. Be- 
sides, the exponent of the density distribution also leads to 
diverging mass with radius and, hence, to the self-gravity 
problem. 

In this paper, we relax the demand of having similarity 
at the ouset and consider instead simple uniform initial con- 
ditions for the gas. In other words, the initial configuration 
is a solid body with a spherically symmetric mass distri- 
bution (a ball) placed in a homogeneous gas with pressure, 
and the problem is to analyze the subsequent evolution of 
this gas under the body's gravity. The evolution will consist 
of infall of gas with a progressive modification of the gas 
distribution around the body, while the gas stays homoge- 
neous far from it. This is essentially the type of accretion 
treated by Bondi (1952), although the artificial imposition 
of stationary flow allowed him to dispense with differential 
equations. At any rate, a homogeneous distribution of gas 
is just a particular case of power law (w = 0) and, there- 
fore, comparison with the results of Chevalier (1989) and 
Kazhdan & Murzina (1994) will be possible. 

These simple initial conditions are suitable for an ana- 
lytic treatment, employing the full power of the methods of 
nonlinear fluid dynamics. These allow us to attain a clear 
picture of the initial non-trivial process near the surface. 
In particular, the exact description of the formation of a 
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shock wave may illustrate its general features in accretion 
processes. Even though the evolution is not self similar, sim- 
ilarity arises in the course of it, and we shall see how and 
why. Due attention is also paid to the self-gravity question, 
and to the decay towards a final static state. 

This paper is organized as follows. In section 2, we in- 
troduce the relevant magnitudes in the problem and, hence, 
various space and time scales, to reach a rough intuitive 
idea of the physics involved and to define the limits of ap- 
plicability of the model. In section 3, wc introduce the fluid 
equations to be used. In section 4, an exact solution is found 
for the initial non-trivial dynamics, which occurs near the 
ball's surface. In the following section wc obtain two similar- 
ity solutions valid for larger t, the first still confined within 
a short distance from the ball, while the second is valid for 
large radius. Finally we consider the long-time asymptotic 
static state and discuss the results. 

A note on notation: wc shall use frequently the asymp- 
totic signs ~ and ~; for example, f{x) ~ g{x) or f{x) ~ g(x) 
(sometimes without making explicit the independent vari- 
able a;). The former means that lim/(a;)/p(a;) when x goes 
to zero or infinity, as the case may be, is finite, while the 
latter moans, in addition, that the limit is one. On the other 
hand, the sign ~ may appear in some instances with the 
related but looser meaning "equal up to a numerical factor 
of the order of unity" . 



2 RELEVANT MAGNITUDES AND 
DIMENSIONAL ANALYSIS 

The initial condition (the solid ball in the homogeneous 
gas) is characterized by four parameters, namely, the ra- 
dius R and mass M of the ball, and the pressure Pq and 
density po of the gas (we assume that the gas is perfect, 
inviscid, and non-heat-conducting or polytropic). In addi- 
tion, we have the constant of gravity G. From these five 
dimensional characteristic paxameters, we can form two in- 
dependent non-dimensional numbers, namely, the ratio of 
densities (4/3)7ri?-^(p„/M) and the ratio RPo/{poGM). The 
quantity Po /po ^ cf) (co being the sound speed) is approx- 
imately the gas thermal energy per unit of mass. On the 
other hand, GM/R is the gas potential energy per unit of 
mass on the ball. Their ratio measures the relative strength 
of gravity in our problem. To have any significant gas infall, 
we must demand that the ball's gravity dominates over the 
gas thermal energy, that is, Rc^/{GM) <t. 1. The ratio of 
densities is also very small. This is a necessary condition for 
neglecting the self-gravity of the gas, as we will do, but it is 
not sufficient. We shall discuss the sufficient condition after 
analysing the typical length and time scales. 

The basic length scale is i?, of course. We can get 
a larger length scale by dividing by the small number 
Rcl/{GM): 



R 



n : 



R 



GM 



Rcl/(GM) 



(1) 



This length TZ marks the scale at which the gas thermal 
energy is similar to its potential energy. It was introduced 
by Bondi (1952) to define a non-dimensional radial variable. 
Alternatively, we may divide R by the cubic root of the other 
small number (4/3)7rfl^(po/M): 



Ji[(4/3)7r(po/M)]i/3 



_3_M 

An po 



1/3 



(2) 



in which we may suppress the numerical factor. Obviously, 
this second length is the radius of a volume of gas such that 
its mass is similar to M. 

Analogously, we have a basic time scale, namely, 
R^^"^ j \fGM , defined only in terms of the ball's parameters 
(and interpreted as the typical time of Keplerian motion 
close to the ball's surface). Dividing by \R^I{GM)fl'^ we 
cancel the R dependence and obtain GM/(^^ that is, the 
time of sound propagation over the distance TZ (note that 
the time of Keplerian motion at distance 7^is Ti^'^/VGM = 
GM/ci as well). Finally, dividing by [R?{po/M)f'^ , we can- 
cel both the R and M dependences and obtain 1 /\/Gpo, the 
typical time of gravitational collapse of the gas (aside from 
the ball, which might act as a seed for the collapse). 

We will assume that the largest typical length (the ra- 
dius of a mass of g£is M) is much larger than the intermediate 
scale TZ; in other words, we consider the mass of gas in the 
volume of radius TZ negligible in comparison with M. Then, 
analogously, the typical time of gas collapse is much larger 
than GM/cq. This is the precise condition for neglecting the 
self-gravity of the gas, understood in an asymptotic sense, 
namely, TZ < (M/po)^/^ or GM/cl < l/^JTTfHi. In non- 
dimensional form, Co/(G'^Af^p()) ^ 1, which can bo restated 
in terms of the two previously defined non-dimensional num- 
bers: besides that both must be very small, the ratio 



R^po ( Rcl \ ^ (TZ\ 
\GM ) \r) 



M 



Tzy^ Wpo 

M 



<1, 



(3) 



that is, the ratio of densities must be much smaller than the 
other one. Interestingly, under this condition, the parame- 
ters M and G will not appear independently in the equa- 
tions of motion but only in the combination GM , so that 
one has only four dimensional parameters. Correspondingly, 
only one non-dimensional number is relevant, namely, the 
ratio R/TZ. 

As regards mass scales, the basic mass is M , of course. 
The mass of gas enclosed in the sphere of radius TZ, namely, 
M. w (4/3)7r72.^po, can be interpreted as the total mass of 
initially bound gas. We observe that M/M ~ TZ^po/M -C 1, 
according to condition (3). On the other hand, this condition 
can also be written as Mj ^ M, introducing the gas Jeans 
mass Mj ~ cl/{G^pof^^. This is, of course, the natural 
mass scale for gas collapse due to self gravity. 

It is convenient to remark that, since we still have after 
neglecting the gas self gravity one non-dimensional number 
(R/TZ), we can construct many length and time scales (large 
or small), but they have no particular physical meaning. 
However, multiplying the basic length 7? by that number, 
we obtain R^Co/{GM) = c^)/ g, where g is the gravity on the 
ball's surface. This small length scale and its associated time 
scale cts/g will play a role in the sequel. 



* With an initial power-law density distribution p = Br 
analogous distance is (M/B)i/(3~"). 
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3 FLUID EQUATIONS 

Under the assumption of spherical symmetry, we are led to 
solving the partial diflterential equations of fluid dynamics 
in one dimension (the radial distance). These equations are 
nonlinear and no general method of solution is available. 

However, given the simplicity of the initial and boundary 
conditions, many results can bo obtained by purely analytic 
means, as we shall see. 

Wc consider an adiabatic evolution of the gas or, more 
generally, a polytropic equation of state, P <x. (7 > 1). 
Then, wc have the contirmity equation, the Euler equation 
and the thermodynamic equation: 



d{pv) ^ 2(H ^ 0^ 



dt dx X 
dv dv 



1 dP 
p dx ' 
d_P_ 

dtp^'^'^d^'p^ ~ ' 



(4) 
(5) 
(6) 



where g = GM/x^ is the gravity acceleration (Chevalier 
1989; Kazhdan & Murzina 1994). The initial conditions are: 
p{Q,x) = po, P{0,x) = Po and v{0,x) = {x > R). The 
boundary condition at the solid surface is v{t, R) — 0. 

In principle, Eq. (6) has the trivial solution P/p'' — 
Po/pQ. Let us introduce the sound velocity c, which satisfies 



c^(p)/co = (p/po) 



7-1 



(C(, = 7P0/P0). 



Hence, the mathematical problem boils down to solving 

dp ^ d{pv) ^ 2{pv) ^ ^ 
dt dx x 



dv ^ ^dv _ 
dt dx 



c\p)dp 



(7) 

(8) 
(9) 



p dx 

The polytropic gas flow in one dimension is amenable 

to powerful mathematical methods (Courant & Friedrichs, 
1948; Landau & Lifshitz 1987; Chorin & Marsden 1993), al- 
though the presence of gravitation complicates the problem. 
However, restricting our interest to a small zone near the 
surface, we can take the acceleration of gravity g constant 
in Eq. (9) (Bisnovatyi-Kogan, Zcl'dovich & Nadezhin 1972). 
This will allow us to take advantage of those methods. 



4 INITIAL STAGES 

Initially, the gas will start falling with acceleration g, that 
is, with a negative velocity increasing as v w —g t. However, 
it will be stopped at the ball's surface, where the density 
and, therefore, the pressure must increase. Consequently, a 
wave transmitting the boundary condition v{t,R) = will 
propagate outwards. It is therefore crucial to determine the 
law of propagation of waves in the present conditions. This 
is called, in mathematical terms, the analysis of character- 
istics. Once this is done, and as long as the dynamics con- 
sists of the propagation of a simple wave, one can obtain 
an exact solution. It is not possible here to describe in de- 
tail how this solution is obtained and our focus will be the 
origin of the shock wave. We refer the reader to Courant 
& Friedrichs (1948), Landau & Lifshitz (1987) or Chorin 
& Marsden (1993) for a general treatment of the theory of 
one-dimensional gas flow. 



For the moment, we confine ourselves to a spherical shell 
over the surface of height much smaller than R, where we 
can consider g constant. Hence, the problem becomes that 
of the fall of gas on a flat surface (the ground) and it is 
convenient to take this surface as the origin of x. So we must 
use the one-dimensional form of the vector divergence and, 
therefore, we must neglect the last term on the left-hand 
side of the continuity equation (8), writing it as 

dp d{pv) 



tt + 



dx 



0. 



(10) 



The gas still unperturbed by the wave merely falls with 
velocity v = —gt. In a reference system falling with it, it 
is at rest and, therefore, the speed of sound is cq. Hence, 
in the original reference system the wave front is located at 
Xf — Cot — gi? 12. Then, for x > Xf, the solution is trivial, 
and wc only need to find the solution for x < Xf. In the 
falling frame, the acceleration of gravity vanishes and the 
set of two equations (10) and (9) reduces to a "gas tube 
problem" that can be solved by introduction of the Riemann 
invariants (Courant & Friedrichs, 1948; Landau & Lifshitz 
1987; Chorin & Marsden 1993). Let us denote x',v',c' the 
coordinate and variables in this frame. Initially, we have a 
constant state, which is preserved for x' > x'f = cot. The 
backward characteristics crossing the forward characteristic 
x' = Cot transmit a constant value of the Riemann invariant 
J_, so the solution is a forward simple wave. Given that 

r I 2c' 2co 

J- =v = -, 

7 — 1 7 — 1 

the sound velocity is related with the gas velocity by c' = 
Co + (7 — l)f '/2- Furthermore, the fact that the solution is 
a forward simple wave implies that the forward character- 
istics arc straight lines. Hence, the wave's propagation law 
v' = F{x — c t) is an implicit equation for v' , assuming that 
we can determine the function F. This is done using the 
boundary condition v{t,0) — 0. The implicit equation is al- 
gebraic and its solution is straightforward. We then obtain: 

-gt, X > Cot — 



V{x,t)={ -E^-l^gt+^(^^ + l^gtY^?^X, (11) 

X < Cot — §f^. 

For a simple wave, the density is a definite function of the 
velocity. We can calculate it from Eq. (7) to be 



p{x,t) = po 



1 + 



7-1 v{x,t)+gt 



Co 



(12) 



These expressions for p and v constitute the solution of 

Eqs. (9) and (10) with the given boundary conditions. 

The alert reader may have noted a peculiarity of the 
solution found above: the wave front's coordinate Xf be- 
gins increasing but, for t > co/g, changes to decreasing and, 
eventually, returns to the origin (free fall with initial veloc- 
ity Co). On the other hand, a simple analysis shows that a 
singularity occurs before, namely, a shock wave. One can 
calculate {dx/dv)t, and its null locus defines a line in the 
a;t-plane, namely, x = [2co -|- (7— l)gt]'^ /(Syg). Initially, the 
corresponding x is larger than Xf, so it is unphysical, but, 
for t = 2co/[(7 + l)g], it meets the wave front; that is, the 
falling gas overtakes the gas just below it and a shock wave 
develops henceforth. From that moment onwards, the so- 
lution (11,12) becomes multivalued and ceases to be valid 
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Figure 1. Evolution of the velocity in the early stages in non- 
dimensional units (co = ff = 1) for 7 = 7/5, displaying how it 
becomes multivalued, hence giving rise to a shock wave. 



(sec Fig. 1). Although it is possible, in principle, to study 
the situation in which a shock wave propagates, the analysis 
cannot be done in terms of a simple wave and it is far more 
complicated. 

Of course, the foregoing analysis is valid as long as ev- 
erything happens in a thin shell, which implies that the typ- 
ical length Co/g <^ R, equivalent to 7?. S> (the already 
mentioned condition of gravity dominated evolution). 



5 SIMILARITY SOLUTIONS 

A one-dimensional gas flow problem may be formulated with 
similarity variables [in which the fluid partial differential 
equations become ordinary differential equations (ODEs)] 
only when the initial and boundary conditions are suffi- 
ciently simple (Sedov 1982). Certainly, this is not our case, 
for we can construct independent length and time scales 
(see section 2). Nevertheless, it is commonly observed that 
nonlinear equations that do not have similarity solutions 
at the outset develop them in an intermediate asymptotic 
regime, that is, a regime between two very different scales 
where both scales can be neglected (Barenblatt 1996). We 
will actually find two independent intermediate asymptotic 
regimes. 

We shall study first the self similar solution that devel- 
ops from the exact solution found in section 4, that is, with 
constant gravity. The corresponding ODE have been formu- 
lated by Bisnovatyi-Kogan, Zel'dovich & Nadezhin (1972). 
However, we shall derive a slightly different (though equivar- 
lent) system of equations more amenable to analytic study. 

Afterwards, we proceed to consider the variation of 
gravity with radius, connecting with the similarity solu- 
tions of Chevalier (1989) and Kazhdan & Murzina (1994). 
These solutions are for an initial density distribution given 
by p = B but a constant density is just a particular 
case. The most interesting solutions is the one with zero 
mass flux near the origin, named "of type 3" (Kazhdan & 
Murzina 1994), since it is expected to match the self-similar 
solution with constant gravity. 



5.1 Constant gravity 

In this case, we have one less parameter, because the pa- 
rameters R and GM only appear in the combination g — 
GM/R?. However, we can still form the two independent 
non-dimensional variables gt/ay and px/co, so we will have 
to remove another parameter to have a similarity solution. If 
we assume that the initial gas temperature is very low, then 
Co — > and the only possible non-dimensional variable is 
^ = x/{gt^) (or a function of it). (We measure the distance 
X from the ground.) This corresponds to the intermediate 
asymptotics Co/g <iC x <C -R; cquivalcritly, in terms of the 
non-dimensional height variable x = gx/c^, it corresponds 
to K i <C n/R (recall that 7^ > Ji). 

To take into account the presence of shock waves and 
the consequent dissipation, we consider the full set of equa- 
tions (5), (6) and (10). We can express them in non- 
dimensional form by introducing new variables: 



X 



p = por, P ■■ 
Some straightforward algebra then yields. 



^Pop. 



(13) 

(14) 

(15) 
(16) 

(17) 



(corresponding to changing to the falling reference frame), 
the term in the second equation vanishes and the system 
of ordinary differential equations simplifies to 



^[u +{u-2)-] + u = 0, 

(w - 2) ^w' -h - w = -r' - r-"' (2p + ^p'), 
(w - 2) ^[\n{py)]' + 2{u - 1) = 0. 

Upon making the change of variables 




+ {u- 2)—] +5 = 
{u — 2) ^-S' + ii^ — u = — r 



(n- 



(2p + 4p), 
2)Clln(p/r^)]' + 2(«-l)=0. 



(18) 

(19) 
(20) 



Bisnovatyi-Kogan, Zel'dovich & Nadezhin (1972) do not 
make the change to the falling reference frame, and restrict 
themselves to solving the equations by a series expansion 
near the ground. Instead, after the change of reference, we 
can follow the general (and more powerful) methods exposed 
by Sedov (1982). In particular, we note that in terms of the 
variable r = In^, the previous equations constitute a system 
of autonomous nonlinear ODE, namely, 

du p{t) [2 -h 7 u{r)] - r{r) u{t) (2-3 u{t) + uirf) 



dr 



dr 
1? 



■ r(r) 



-yp{T)+riT)[u{T)-2Y 

-2p(r) -I- r(r) [w(t) - 2] w(t) 

(-7P(r) + r(r) [{i(r) - 2]^) [u(r) - 2] 



(21) 



, (22) 



dp ,^ 27p(r)- r(r) (4- (6 + 7) '"(r) + 2w(r)') 

— = p{t) ^ 5 -. (23) 

dr ^ _^p(^)+r(T)[M(T)-2]^ ^ ' 

One can study this system of equations with the methods 
of nonlinear ODE, that is, study its singular points, phase 
portrait, etc. However, given the homogeneity properties of 
these equations with respect to p and r, it proves convenient 
to introduce the variable 9 = p/r. It is (in the falling frame) 
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the non-dimensional form of the temperature for the perfect 
gas with pressure P and density p, that is, T = fj, (x'^/t^) 
{jj, being the molar mass and T being measured in energy 
units). Hence, 



de^^^e (1 + 7(m 
du 



■2)) - (m-2) (4- (5 + 7) m + 2^2) 



(m-2) (fii (2 + 7m) -m (2-3m + m2)) 



(24) 



The solution of this equation provides the relation between 

the "velocity" u and the "temperature" 6, 9(u). Substituting 
it back into Eq. (21), we have an ODE for m(t), which is 
immediately solved by a quadrature. Analogously, one solves 
for r(T). 

To solve the nonlinear ODE (24) , let us note that we are 
actually interested in some particular initial conditions, do- 
rived from the initial and boundary conditions of the original 
partial differential equations (PDE). The initial conditions 
for the PDE are given at t = 0, that is, at ^ = ^ = oo. They 
are p(0, x) = po, P{0, x) = Po = (co — > 0) and v{0, x) = 0, 
implying that r = 1 and u = 0. The boundary condition for 
the PDE is given at a; = => ^ = and f = 1/2, where 
V = 0. Since v = gt{£,u — 1), it implies that u = 2. 

Let us see if we have gathered enough information to 
determine a unique solution of Eq. (24) . Given that ^ — > as 
^ ^ oo, this must be a fixed point. There are three singular 
points of Eq. (24) with ^ = 0, namely, with u = 0, 1 or 2, 
respectively. If we want to have a finite w as a; — > oo, we 
must take the point (0,0). So we know that the solution of 
Eq. (24) that we seek must end at the origin. On the other 
hand, the singular point (2, 0) corresponds to the boundary 
conditions at a: = 0. However, the solution that departs from 
this point does not end at the origin (see Fig. 2). This does 
not mean that there is an inconsistency, for there can be 
discoritirmities in the solution. Indeed, in most self-similar 
solutions the boundary conditions can only be fulfilled if 
there is a point of discontinuity (Sedov 1982). 

Therefore, recalling the exact solution of section 4, 
wo observe that a shock wave arises and propagates up- 
wards. Then we must consider discontinuities associated to 
the presence of a shock wave. The conditions for conserva- 
tion of mass, momentum and energy across a shock surface 
read (Courant & Priedrichs, 1948; Landau & Lifshitz 1987; 
Chorin & Marsden 1993) 



plVl = p2V2, 

Pi + piv'i = P2 + p2vl, 



(25) 



-I- Wl 



+ W2, 



where subindices refer to the values on each side of the shock 

and w is the enthalpy (or the appropriate thermodynamic 
function for general 7); for a perfect gas, w = 7P/[p(7 — 1)]. 
These equations hold in a coordinate system in which the 
shock surface is at rest. On the other hand, the location of 
the shock must be at fixed ^, say ^s- Consequently, the shock 
wave velocity is 



dxs dt r»/- . 0*^3 



that is. 



= 2. Then, 



ri(Mi - 2) = r2{u2 - 2), 

Pi + ri{ui - 2f =P2 + r2{u2 - 2)^ 



(ui-2y 



(26) 



, 7 PI 

T— 1 ri 



_|_ 7 P2 



0.4 



0.2 




0.5 



1.5 



Figure 2. Relevant solution of the ODE (24) for 7 = 7/5, dis- 
playing the origin and its image at the shock surfa<;e, the point 
(5/3,5/9). 



valid in both the rest and the falling reference frames. Using 
the falling frame and eliminating r2 /ri between the first and 
second equations. 



{ui — 



+ ui-2 = 
2f + ^ 



92 



?1 = (U2 



+ U2 - 2, 

2 I 2t_ 



2Y + ^fe. 



(27) 



We have arrived at an involutive mapping of the plane 
(u, 9) as the relation between the variables at either side of 
the shock discontinuity. Since we have the condition that as 
a; — > 00 the solution goes to (0, 0) and (in the falling frame) 
these values hold all the way down to the shock surface, we 
just need the image of the origin under the mapping. This 
yields the point 



4 8 (7-1) , 
7 + 1' (7+1)' ' 



Therefore, we need the solution of Eq. (24) that goes through 
this point. We can see in an example (Fig. 2) that the strand 
of this solution that we need ends at the point (2,0), hence 

satisfying the boundary condition at a; = 0.^ 

Now, substituting the function just computed back into 
Eq. (21), we obtain 



dr 



-i0{u) - {u-2Y 



-e(u) {2 + ^u)+u (2-3u-|-u2) 



du. 



(28) 



from which we can deduce the interval of r between the 
points (4/(7 -I- 1), 8 (7 - 1) /(7 + 1)^) and (2,0) by integra- 
tion. The result for 7 = 7/5 (the adiabatic index of a per- 
fect diatomic gas) is 0.0880104. Hence, we obtain the quo- 
tient between the corresponding values of 5. This quotient 
gives the location of the shock relative to the ground: re- 
covering the original notation that distinguishes both refer- 
ence frames, fs/(l/2) = exp 0.0880104 = 1.09200 ^ = 
0.0919995/2 = 0.0459997, so that the coordinate of the 
shock is Xs = 0. 0459997. c/f^ 

It is easy to compute the post-shock u in the rest frame: 



^ Notice that this implies that the asymptotics co — > is of 
the first kind (the simple case) according to the denomination 
of Barenblatt (1996). 
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5.2 Variable gravity 




-100 V 



In X 



Figure 3. Self-similar velocity plot (cq = g = 1, 7 = 7/5) show- 
ing two sets of lilies, for constant t or ^, respectively (the border 
corresponds to the shock at ^a) . The v profile at constant t is 
transported homologously along the constant ^ lines. 



u = [(1.092/2) (5/3) - l)]/0.046 = -1.957. Hence, the post- 
shock velocity is v = {xs/t)u = ^sugt = —0.095*. Com- 
paring with the pre-shock value v = —gt, we see a great 
reduction. The corresponding kinetic energy is dissipated as 
heat. The similarity properties of the velocity and its jump 
at the shock are shown in Fig. 3. 

The ratio of post-shock to pro-shock densities can bo 
derived from the first Eq. (26) as pilpi = (7 — l)/(7 + 
1) = 6, as corresponds to "strong shocks" (Landau & Lif- 
shitz 1987). The post-shock temperature is given by T = 



p(xlit^)e = p (8 (7 - 1) /(7 + 1)') lla't^ 



0.166 /lig^r 



T/To = 0.232 {gt/caf. Therefore, the post-shock pressure, 
P = pT/ IX, is not polytropic. 

Other properties of the self-similar solution found are 
worth noting, but require further analysis of the ODE sys- 
tem (21, 22, 23). For example, at the ground, 6 = C^^/? = 
0, so that 9 = limj^o S/i^^^)- To calculate this limit we need 
to know the behaviour of near ^ = (r = — In 2). The ex- 
pansion of the system (21, 22, 23) yields 6'^'^^ ~ ■& — 2 w 
-2(1 + l/7)(r + In 2), r ~ ^"^ with S = 1/(7 + 1). This 
implies that w(0) = -2/7, 6 Furthermore, p ~ 

which means that P = {x/t)'^ppo oc g^t^po at a:: = 0. The 
density diverges as a: — > but is integrable (as is necessary 
to have a finite mass). Since the self-similar solution is not 
valid if X < c^/g, wc can interpret this divergence as result- 
ing from great concentration of gas in the region x ~ Co/g, 
relative to its initial value. In fact, it is easily derived that 
this mass ratio grows as (gt/co)^^ ■ 



If we consider heights of the order of the planet radius or 
larger, we have to return to the full continuity equation (4) 
and reset g = GM/x^ in the Euler equation (5) (a: is again 
the distance from the center). We can find similarity solu- 
tions in the intormodiato asymptotics R <^ x <^ TZ, in terms 
of the non-dimensional variable ^ = x^ / {GMt^). Hence, we 
can express the continuity and Euler equations as 



(3m - 2) ^w' + u 



C[3w'-f (3m-2)^]-F3m = 0, 
-M = -r'-r-'(2p + 3^p'), 



The energy equation becomes 

(3u - 2) ^[ln(p/rT)]' + 2{u - 1) = 0. 



(29) 
(30) 

(31) 



These three equations look similar to the ones cor- 
responding to constant gravity, and they can also be re- 
duced to two differential equations for ii(^) and 6(^), re- 
spectively (Cheng 1977, Kazhdan & Murzina 1994). Unfor- 
tunately, they are more difficult to analyze: the change of 
variables that transformed the constant-gravity equations 
into an autonomous system is not available. In fact, the free- 
fall (pressureless) problem is now nontrivial, although it can 
be solved. 

If p = 0, Eq. (30) decouples, becoming an equation for 
just m(C), namely, 



(3'u — 2) + u' 



(32) 



To implement the PDE initial condition, at ^ —» 00, it is 

convenient to make the change of variable ^ = trans- 
forming the equation into 



(33) 



du _ u — w -|- ^ 
di ~ (3« - 2) I ' 

The point (0,0) is a singular point of the ODE: a solution 
is unspecified, unless we introduce am additio nal con dition. 
The initial condition v = implies that v = y/ GM/x ^"^^^u 
goes to zero as ^ ^ 0. To impose it, it is best to linearize 
and solve the ODE around (0,0), obtaining 



(34) 



-1/2, 



so that the singular point is a node. Since lim|_^(j ^ 
C, the particular solution of the nonlinear equation (33) in 
which we are interested is the only one that satisfies u'{0) = 

-1 (C = 0). 

The full free-fall solution can be obtained in parametric 
form with the Lagrangian formalism (see Appendix 1). It 
reads 



8 cos a 



[2a + sin(2a)]2 ' 



, sin(2a) , sin a 



9 cos a — cos(3 a) -|- 12 a sin a ' 



(35) 

(36) 
(37) 



where a G (0,7r/2). The corresponding graphs for w(^) and 
r(^) are shown in Fig. 4. 

This free-fall solution is to be distinguished from the 
simpler free-fall solution considered by Cheng (1977) or 
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Figure 4. Similarity solution for free fall with vari- 

able gravity. 



Kaahdan & Murzina (1994), namely, u = — y2/^, which 
corresponds to vanishing initial energy or, in other words, 
to the initial configuration consisting of gas particles at rest 
at infinite distance. Note that t he d ominant term of the so- 
lution (34) is proportional to and coincides with their 
free-fall solution for C = —\/2, being then an exact solu- 
tion of the nonlinear equation (33). Furthermore, for other 
non-null values of C, we have different self-similar asymp- 
totic solutions, with a given initial velocity proportional to 
^JGmJx. On the other hand, we can derive u = —\f2i\ 
from the parametric solution in the limit C ~* (correspond- 
ing to Q ^ 7r/2 and meaning small distance or long time), 
since the initial energy of the gas particles becomes negligi- 
ble in this limit. 

The complete similarity solution consists of the free-fall 
solution and an inner solution with pressure matching across 
a spherical shock wave. The PDE's boundary conditions at 
a: = determine the form of the inner solution of the two 
diflferential equations for m(^) and ^(^) near ^ = 0. It must be 
a solution of "type three" , with one free parameter (Kazhdan 
& Murzina 1994): 



where uo and are related by 



(38) 



[3(7 - l)/3 -I- 37 - 4]wo = 2/36io 



and /3 = -1/3-^(2/3)^(7- l)/7. The match with the pre- 
vious free-fall solution across the shock determines the free 
parameter. There is no straightforward way to find the shock 
coordinate and a sort of "trial and error" procedure is 
necessary (Chevalier 1989, Kazhdan & Murzina 1994). An 



approximate value can be obtained by matching the expres- 
sions (38), namely, = 0.013. 

Both the pre-shock and the post-shock gas velocities are 
given by v{;xs,t) = (GM^s)^^^m(6) t''^''^, for the respective 
u values. Hence, the similarity mass flow across the shock 
is m = 'inx^pv = 47rGMpoCs'"(^s)^(^s) which increases 
with time. This may seem paradoxical, since the velocity 
decreases as t~^^^ , but the increase in shell mass due to the 
spherical geometry compensates for it. 

We note that the shock velocity experiences at a:^ ~ Ji 
a crossover from Vs = 2xs/t = 2^sgt to Vs = 2/3 (xs/t) = 
2/3{GM^sy/^t-^/^ (with different ^s), turning from ac- 
celerating to decelerating. The previous solution is valid 
for Xs — R <^ R, whereas the present solution is valid for 
Xs ^ R. Once the second solution takes hold, given that 
the sound velocity for the far-away gas is co, the shock dis- 
appears after its velocity becomes subsonic. This occurs for 
{GM£_sf/^ t-^/^ CO ^ t GM/4 and Xs ~ TZ. 



6 THE STATIC STATE 

Since the shock wave eventually vanishes, leaving the gas 
with more entropy than in the initial state, we expect that 
it must roach a stationary state in the long run. The con- 
dition that v{x) = (from the continuity equation and the 
boundary condition) implies that the stationary solution is 
actually static. 

The hydrostatic equation, 

— = ~ gdx, 
P J 

with P cc p'' and g constant, has the solution: 

7-1 



(39) 



where the subscript b denotes values at the ball's surface. 

Actually, this formula describes well the lower region of the 
present Earth's atmosphere. Nevertheless, it is unphysical 
when taken over the entire range of x (for 7 > 1, that is, for 
a non-isothermal distribution): p'^~^ (which is proportional 
to the temperature) decreases linearly with height, becoming 
negative for some height. It is natural that there cannot be 
a stationary state with constant gravity, given that nothing 
can prevent the free fall of gas. 

In contrast, if we take variable g = GM/x^ , 



^ =1- 



7 - 1 pbgbR 
7 Pb 



('-!)■ 



(40) 



In this case, it is more convenient to take the reference at 
infinity, where the density and pressure keep their initial 

values, so that 



T_ 

To 



7-1 



1 + 



7 - 1 po GM 
■y Po x 



l + (7-l)|. (41) 



As expected, TZ gives the sphere out of which the gas can 
be considered gravitationally unbound (for its thermal plus 
gravitational energy is positive) and, consistently, the scale 
of crossover to the asymptotic {x 00) density po. 

For small radius a: <C 72., we can neglect the unity in the 
right-hand side of Eq. (41), so that we have power laws. In 
particular, the density can be written as 
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^=(^r, (42) 

Po \nxJ 

by introducing n according to 7 = 1 + 1/n. It coin- 
cides with Cheng's hydrostatic solution, after substituting 
poTT- = (GM)"7-"pS'"^iPo-". However, Cheng's subscript 
indicates arbitrary reference values, given that p"+ip~" is 
constant. Therefore, we see that it is convenient to take the 
asymptotic values as reference, as long as they are non van- 
ishing, Ti, being the scale where p and P take those values. 
Moreover, we can let those asymptotic values vanish, by tak- 
ing TZ ^ 00 while the product p^TZ" has finite limit. On 
the other hand, since the ratio between the surface values 
of temperature or density and the initial values is given by 
Tb/To = (pt/po)^^" = TZ/{nR) ^ 1, we may prefer to choose 
the former values as reference. Then we can analogously take 
the limit R and pt 00, while keeping pbR" constant. 
Since we have pure power laws when 71/ R 00, the refer- 
ence is arbitrary. 

We must remark that the law (41) gives a density con- 
trast p — Po that decreases too slowly with the radius, yield- 
ing a divergent accumulated mass of gas as a; ^ cx). In fact, 
the exact form of the hydrostatic equation involves the grav- 
ity of the gas and, hence, is equivalent to the equation that 
describes stellar structure, namely, the Lane-Emden equa- 
tion (Chandrasekhar 1967). If the bound mass of gas is 
much smaller than the mass of solid material, that is, if 
TZ <ti {M / poY^^ , the solution of this equation is well ap- 
proximated by Eq. (41), while for x > {M/poY^^ the decay 
is faster than the decay ~ 1/x given by it. One can then 
approximate the mass of accreted gas by the total mass of 
gas inside the sphere with radius TZ, which is larger than the 
initially bound gas mgiss, M., because a portion of the gas 
that is initially out of the sphere falls inside. The detailed 
calculation in Appendix 2 shows that the ratio of this por- 
tion to M is of the order of unity for 7 < 4/3 (in particular, 
it takes the value 2.65 for 7 = 7/5)- 

We have seen in the previous section that the spherical 
shock vanishes when it reaches a radius ~ TZ. This is con- 
sistent with the fact that the gas outside a sphere of radius 
TZ is inmune to the gravitational influence of the planet. On 
the other hand, since the shock vanishes for t ~ GM/cq, 
this implies that the similarity solution becomes useless for 
t > GM/cq, marking the crossover to the static state. Nev- 
ertheless, it is necessary to remark that the static state is 
an intermediate asymptotic regime not only in space, as was 
commented above, but also in time, being always required 
that t < l/s/GpQ. 



7 DISCUSSION 

We have performed a thorough analysis of a simple model 

for gas accretion onto a solid body. We have seen that there 
are two relevant length scales, namely, TZ and {M / poY^^ , 
the relative value of which determines whether the dynamics 
consists of accretion or collapse: to neglect self gravity, it is 
required that TZ < {M/pof^^. 

We have not restricted ourselves to similarity solutions. 
We have divided the space dimension (the radius x) into a 
near zone (near the ball's surface, where x — R <^ R and g 
is constant) and a far zone {x ^ R). In the near zone, we 



have obtained an exact solution for short time and a self- 
similar intermediate-asymptotic regime for longer time and 
radius. In the fax zone, we have a self-similar intermediate- 
asymptotic regime for relatively long time and radius and an 
exact static state for very long times. These solutions pro- 
vide a good intuitive understanding of the entire dynamics. 

The crucial feature is the formation and propagation 
of a shock wave until its dissipation. The shock wave arises 
in a short time ~ co/g and it rapidly grows. As long as 
its propagation is strongly supersonic, we can consider the 
gas cold (co = 0) or, in other words, the dynamics driven 
by gravity. The shock wave experiences a crossover between 
one stage of strong dissipation and acceleration to another 
of slowing down, owing to decreasing gravitational energy 
input, until its eventual vanishing, leaving behind dense and 
hot accreted gas in a quasi-static distribution. 

The process of accretion begins at t = and contin- 
ues until the end, when v ^ 0, but it can be considered 
finished when t ^ GM/cq. Since the velocity experiences a 
strong reduction at the shock, with a consequent increase in 
the density, it is sensible to define the accretion rate as the 
value of the mass flow rh there. We have seen that this value 
is proportional to t during the second self-similar regime, so 
the amount of mass being accreted only begins to decrease 
at the end of it. The total accreted mass in the static state is 
not very large if 7 > 4/3 (7 = 4/3 for the black-body opac- 
ity limit), in spite that the corresponding ratio of surface 
density to po is very large ([TZ/(nR)]^). 

Although an adequate polytropic index accounts for 
some form of heat conduction or radiation, the polytropic 
static state for a heat conducting or radiating gas is unsta- 
ble against further cooling. Therefore, we must remark that 
the true static state, from a strictly theoretical thermody- 
namic point of view, is only achieved by cooling until reach- 
ing an isothermal state. The corresponding isothermal den- 
sity distribution is much more concentrated near the ball's 
surface (it is exponential), leading to a far larger accreted 
mass. However, the time to reach this isothermal distribu- 
tion would be, arguably, larger than the gas collapse time. 

We may wonder how our picture would be modifled by 
different initial conditions, namely, by a non-homogeneous 
density distribution (e.g., a power law). Gravity only hcis 
effect inside the sphere of radius TZ, in which the gas rapidly 
becomes supersonic and approaches free fall, until the shock 
wave arrives. The free-fall solution (35,36,37) is independent 
of the value of po, which can be a function of x (see Appendix 
1), but the form of the shock wave and the post-shock gas 
distribution will change. 

We may also consider an initial non-homogeneous tem- 
perature distribution. Near the body surface, we can make a 
linear approximation of it. If its slope is sufficiently smaller 
(in absolute value) than the static value provided by Eq. 
(39), we can still consider it isothermal and predict that the 
fall of gas will produce a shock wave. Nevertheless, we can 
envisage distributions close to the static one such that the 
gas falls gently on the surface without ever giving rise to 
a shock wave (such as, e.g., the distribution corresponding 
to the final stage of our process, after the vanishing of the 
shock wave). Presumably, these distributions lead to little 
accretion. 

Now, I briefly estimate the numerical scales suitable for 
the application to neutron star accretion or planet accretion. 
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For 

-2 



I neutron star of 1 Mq with R ~ 10 m, g ~ 10 m 
s The sound speed (proportional to the square root of 
the temperature) can be large, say co — 10000 m s"\ Then, 



co/g ^ 10- 



s and C(,/g '. 

-,12 



10" 



/GM = y^R/g ~ 

10"* s, but 7^ ~ 10'^ m and GM/cq ~ 10® s. It is not clear 
how to determine po. A "typical" galactic value would give a 
collapse time similar to star formation times, that is, of the 

order of a milhon years. Considering a planet placed in the 
protoplanctary nebula, wo may take for g the Earth value of 
10 m s^"^ and, for example, co ~ 300 m s~^ (its current value 
in the low Earth's atmosphere); the initial typical scales of 
the problem, namely, the time co/g and the distance Co/g, 
take value s of 30 s and 9 Km, respectively. The basic time 
scale is R/g ~ 800 s (for an earth-like planet with R ~ 
6000 Km). We further estimate 7^ ~ GAf/cg ~ 5 10*' Km 
and GM/cq ~ 5 10^ s ~ 1 yr. A plausible value for the 

g/cm^ 



density of the protoplanetary nebula is po ~ 10" 



We see that its collapse time 1/^ 



s is sufficiently 



larger than GM/cq for neglecting self gravity. 



It is clear now that 
r'^^ = F{x/xo), 
where 

F{y) = 



(7) 



-3/2 



i: 



da 



V2(l/^-l)' 



Inverting (7) we obtain x/xo as a function of ^, which upon 
substitution in (4) yields its self-similar form. This is the 
solution of the ODE that satisfies the appropriate initial 
condition (see section 5.2). 

The method exposed above works in general for free 
self-similar motion in an external gravitational field. In the 
present case, one can calculate the integral in Eq. (6) by the 
change of variables a = cos^ (p. This allows one to obtain the 
solution of the diflferential equation (1) in parametric form: 
Let 



cos^ a •■ 



■ x/xo; 



(8) 
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then, 

VGMt 1 , sin(2a), 
and 

c-i/2 1 -3/ , sin(2a) 
(a)(a+ 2 — ^' 

-3/ \/ , sin(2a), . 
w = — cos (a) (a H J- — ^) sma, 



(9) 

(10) 
(11) 



Appendix 1: Lagrangian similarity solution for the 
free fall in the field of a point-like mass 

We describe here the "Lagrangian solution" of Eq. (32), 

(3u-2)5'u' + - It = (1) 

In the Lagrangian picture, the velocity is given by the energy 
equation 



GM 



2 X 



(because the energy is conserved), that is, 



X 



(2) 



(3) 



The initial condition that the particle is at rest implies that 
E = -GM/xo. Hence, 



t 

u = —V 

x 



2C-Hi-f-), 

Xo 



(4) 



and to have a self-similar form we must express x/xo as a 
function of ^. In order to do it, we must solve the equation 
of motion (3): 



t = 



ds 



^2{GM/s - GM/xo) 



(5) 



First, we write it in a self-similar form, by introducing a ■ 

s/xo, 



t ■■ 



3/2 

/GM 



X/XQ 



da 



V2(l/^-l) 



(6) 



where a G (0, 7r/2). It is easy to verify that this parametric 
solution fulfills the initial condition lim{_Kx> u — —1. 

In the Lagrangian formulation, the density is given by 
(Landau & Lifshitz 1987) 



/ dxp \ Xo 
V dx Jt x^ 



po, 



where Xq/x^ is a spherical geometry factor. From Eqs. (8), 
(9) and (10), 



Po 9 cos(a) — cos(3 a) -|- 12 a sin(a) ' 



(12) 



which together with Eq. (10) constitute the parametric equa- 
tions of the density. 



Appendix 2: Accreted mass in the static solution 

Let us calculate the gas mass enclosed in the sphere of radius 
TZ in the final static state: 



■^f = / P{^) 47ra; dx = 

J R 

47rpo / (i + (t-1)7)'' ^^^^ (13) 

This integral is convergent for J? — > if 7 > 4/3. We can 
easily express it in non-dimensional form: 



Mf = a-kpq-r' f (i + i-r- 1)-) ' 



1\t-1 2, 

s ds. 



(14) 
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Figure 1. Accreted mass ratio A as a function of n (7 = 1 + 1/n). 



Now, it is convenient to make the change 7 = 1 + 1/n. The 
resulting integral can be computed in terms of a hypergeo- 
metric function: 



V nsj 



s^ds 



R/7L 

1 



(3 - 7i)n" 



2-Fi (3 — n, — n, 4 — n; —ns) 



(15) 

R/TZ 



We are interested in the limit of [Mf — M]/M when 
R/TZ ^ Q, in which we are left with the value of the bracket 
at s = 1. Denoting A = limi{/7j_>o A^//A^ — 1, 



A = 3 / {1 + ^^"s^ds- 1 = 



3 

2Fi(3-n, -n,4-n;-n) - 1. (16) 



(3 — n)n 

The value of A grows almost linearly, with a small slope, 
up to near the pole n = 3 (7 = 4/3), as Fig. 1 shows. In 
particular, it takes the value 2.65099 for 7 = 7/5 (n = 2.5). 

The polo shows the divergence of the integral as s — »■ 0. 
Then, it is not possible to make R/IZ = near the pole. The 
n = 3 integral yields 

Its value is smaller than 5 even for TZ/R as large as 10^. 
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